Climate change will lead to range shifts and genetic diversity losses of dung beetles in the Gobi Desert and Mongolian Steppe

Desertification is known to be a major threat to biodiversity, yet our understanding of the consequent decline in biodiversity remains insufficient. Here, we predicted climate change-induced range shifts and genetic diversity losses in three model dung beetles: Colobopterus erraticus, Cheironitis eumenes, and Gymnopleurus mopsus, distributed across the Gobi Desert and Mongolian Steppe, areas known for desertification. Phylogeographic analyses of mitochondrial COI sequences and species distribution modeling, based on extensive field investigations spanning 14 years, were performed. Species confined to a single biome were predicted to contract and shift their distribution in response to climate change, whereas widespread species was predicted to expand even if affected by range shifts. We indicated that all species are expected to experience significant haplotype losses, yet the presence of high singleton frequencies and low genetic divergence across geographic configurations and lineages mitigate loss of genetic diversity. Notably, Cheironitis eumenes, a desert species with low genetic diversity, appears to be the most vulnerable to climate change due to the extensive degradation in the Gobi Desert. This is the first study to predict the response of insects to desertification in the Gobi Desert. Our findings highlight that dung beetles in the Gobi Desert and Mongolian Steppe might experience high rates of occupancy turnover and genetic loss, which could reshuffle the species composition.

that could arise in response to climate change 15,16 .However, SDM approaches, which are based on limited environmental data, have been criticized as oversimplistic since they do not account for evolutionary processes and biological features that affect species' distribution 17,18 .For instance, integrating SDMs and phylogeographic information would provide a powerful framework for inferring the effects of future climate change on target organisms 5,19 .
Phylogeography may provide insight into the contemporary and future responses of species to environmental changes.Indeed, genetic structure can provide insight into a species' evolutionary history and into the dynamics of contemporary gene flow among populations, and genetic diversity can be used to evaluate species' evolutionary potential and resilience to environmental change 4,5 .In addition, populations that suffer from unsuitable habitats or limited connectivity often harbor low genetic diversity as a result of bottlenecks or inbreeding depression and are, as a result, highly vulnerable to environmental change 20 .Therefore, investigating the genetic diversity and structure of populations across species' distribution ranges is critical to evaluating the vulnerability and adaptive capacity of species to environmental change 5,21 .
Dung beetles are key organisms in terrestrial ecosystems because they contribute to various ecosystem processes, including nutrient cycling, parasite suppression, and secondary seed dispersal 22 .The group has been widely used as a bioindicator of biodiversity and ecosystem health because of their ecological importance, sensitivity to environmental factors, and established standardized sampling methods 23,24 .Serious concerns exist regarding the global decline in dung beetle diversity, and empirical data show that populations of many dung beetle species have significantly declined in response to habitat fragmentation 25 and climate change 26,27 .
Mongolia spans a vast plateau with wide steppe and desert regions.More than 67 million livestock, including horses, cattle, sheep, goats, camels, and yaks, are farmed across the Gobi Desert and Mongolian Steppe, mostly through nomadic herding (National Statistics Office of Mongolia; https:// www.en.nso.mn/, accessed on April 08, 2022).This environment, abundant in food resources throughout the country, dung beetles are essential for maintaining the ecosystem health.To date, 78 dung beetle species have been documented in Mongolia 28 .Research on the conservation of dung beetles in Mongolia has been primarily confined to a few national reserves 29,30 .However, Mongolia is a particularly well-suited region for investigating the effects of desertification through climate change on organisms.Encompassing two distinct climate regions and less affected by other anthropogenic disturbances such as habitat fragmentation, and chemical treatments, the country act as a natural laboratory to pinpoint and understand the direct implications of climate change.
Accordingly, this study explored how future climate change-induced range shifts and the loss of genetic diversity may differ among dung beetles in desert and sub-arid steppe region by incorporating SDMs and phylogeographic analyses.To achieve this, we investigated the genetic diversity, population structure, demographic history, and distribution of three model species from the Gobi Desert and Mongolian Steppe in Mongolia.We then predicted their future distribution and genetic diversity under climate change scenarios.

Population structure and demographic history analyses
Populations of C. erraticus and G. mopsus exhibited significant pairwise genetic differentiation, whereas those of C. eumenes did not.The pairwise differences (F ST ) of C. erraticus populations ranged from − 0.066 to 0.367 (Supplementary Information, Table S1).However, the genetic differentiation between CER1808 and the other three sites was significant, which indicate that the local population is genetically isolated from the others.Meanwhile, the pairwise differences (F ST ) of C. eumenes ranged from − 0.087 to 0.081 (Table S2), with all F ST values indicating low levels of genetic differentiation.Finally, the pairwise differences (F ST ) of G. mopsus populations ranged  S3).As in C. erraticus, most of the F ST values were relatively small, indicating a low inter-population differentiation level.Only four of the 45 comparisons among the 10 populations indicated significant differentiation.
In C. erraticus, AMOVA indicated significant spatial genetic differentiation between the Khangai Mountains group and other populations (Table 2).In G. mopsus, AMOVA indicated that significant spatial genetic structuring occurred under climate and geography, as reported by previous study 31 .However, no spatial genetic structuring was observed in the geographical groupings of C. eumenes.

Predicting distribution and genetic diversity under climate change scenarios
Estimations of future C. erraticus distribution indicated a westward shift.The northeastern area of the species' distribution contracted and shifted among the various time points and climate scenarios (Fig. 2).Under all climate scenarios, the distribution of C. erraticus exhibited a net reduction in area, ranging from − 18.55 to − 46.26% and from − 31.80 to − 63.23% in 2050 and 2070, respectively (Table 3).In contrast, the future projections of C. eumenes revealed a severe reduction in the southern desert region and a northward shift to central steppe region (Fig. 3).The distribution of C. eumenes was projected to either increase or decrease, with predicted changes in area ranging from − 79.01 to 19.98% and from − 28.02 to 58.42% in 2050 and 2070, respectively (Table 3).For G. mopsus, the species' distribution was predicted to shift northward, with a net expansion in the northern steppe region (Fig. 4).Under all climate scenarios, the distribution of G. mopsus exhibited a net expansion in the area, ranging from 27.23 to 33.49% and from 21.63 to 43.44% in 2050 and 2070, respectively (Table 3).The predicted loss of future genetic diversity in the overall population varied based on the climatic niche and phylogeography of each species.For example, six of the 12 C. erraticus populations were predicted to be affected by climate change (Table S4).The CER1402, CER16EX2, and CER1822 populations were predicted to be extirpated in every scenario, whereas CER1808 (currently exhibiting low suitability) increased in suitability under most of the future climatic scenarios (all but SSP585 in 2070).In addition, 24 (54.6%) of the 44 C. erraticus    www.nature.com/scientificreports/haplotypes were estimated to be affected by at least one climate change scenario, and six were predicted to be completely lost in response to climate change (Fig. S6).For C. eumenes, all populations were predicted to be affected by climate change (Table S4).The CEU22A4 and CEU22A8 populations were predicted to be extirpated in every scenario, and under the most severe scenario, local extinction was predicted for all populations by 2070.
In addition, five (19.2%) of the C. eumenes 26 haplotypes were predicted to be lost in response to climate change (Fig. S7).For G. mopsus, six of the 10 populations were predicted to be affected by climate change (Table S4).The GM1510 population was predicted to be extirpated in every scenario, whereas GM1403 (currently exhibiting low suitability) was predicted to increase in suitability.In addition, 107 (79.3%) of the 135 G. mopsus haplotypes were predicted to be affected by at least one climate change scenario, and two haplotypes were predicted to be lost (Fig. S8).However, despite the predicted losses in the number of haplotypes, none of the three species were expected to undergo large reductions in haplotype or nucleotide diversity (Table 3).Similarly, there was no evident loss of phylogenetic diversity (e.g., loss of clades) in the phylogenetic tree of each species.Instead, the missing haplotypes were scattered throughout the tree.

Discussion
This study is the first to model the effects of climate change on insects from the Gobi Desert and Mongolian Steppe region.Given the serious concerns about desertification in the Gobi Desert and adjacent regions 7,8 , there have been few studies investigating its influences on organisms there [11][12][13] .This issue fundamentally stems from the lack of accumulated data in this region.Our results on the distribution and genetic diversity of Mongolian dung beetles can be valuable data for future studies on global-scale desertification as well as for the responses of organisms to desertification in the Gobi Desert.Genetic diversity and structure can be used as indicators of the tolerance and adaptive potential of a population to environmental stress 32 .In this study, C. erraticus and G. mopsus exhibited high genetic diversity with high singleton frequency, which could reflect the species' adaptive potential.In contrast, the desert-dwelling C. eumenes exhibited moderate genetic diversity, with a single haplotype accounting for 65% of specimens, which could indicate that the adaptive potential of this species is lower than that of either C. erraticus or G. mopsus.
Most of the F ST values among the populations of the three species were relatively small, suggesting low levels of interpopulation genetic differentiation, which could be related to contemporary gene flow or demographic events such as bottlenecks, recent divergence, or selection.For example, there were no significant F ST values for distant populations, C. erraticus populations at CER1822 and CER1502 (over 1100 km), C. eumenes populations at CEU1609 and CEU1713 (over 450 km), and G. mopsus populations at GM1417 and GM2207 (over 790 km; Tables S1-S3).Notably, significant differences were observed between the Khangai Mountains and other populations of C. erraticus (Table 2).These differences could be derived from historical demographic events involving geographic barriers, for instance, the founder effect or selection on historically small, isolated populations from genetically diverse population in the past 33 .This phenomenon was observed only in the Khangai Mountains populations, whereas the Altai Mountains population (CER1822) did not exhibit significant differences from the other populations (except for CER1808).The genetic structures in both geographic and climatic group for G. mopsus coincided with those reported by previous study 31 .
The demographic histories of the three dung beetle species were observed to differ.The haplotype networks of all three species showed major haplotypes surrounded by many unique haplotypes.For C. eumenes, a single major haplotype was identified, while multiple major haplotypes were observed in the other two species.Our demographic history analyses indicated that C. eumenes and G. mopsus have undergone recent rapid expansion.A previous study proposed the domestication-driven expansion hypothesis to explain this recent expansion of Mongolian dung beetles 31 .However, Tajima's D analysis and the mismatched distribution of the overall population revealed that C. erraticus maintained a constant-sized population rather than undergoing a recent expansion.One plausible explanation for this result could be the difference in food ranges among the three species.C. eumenes and G. mopsus prefer the feces of large herbivores, primarily livestock, whereas C. erraticus has a broader food range that includes the feces of small mammals.The presence of C. erraticus inside the burrows of rodents has been reported, indicating that the species can utilize the dung or residues found within burrows 34 .Therefore, C. erraticus may have been less affected by domestication events than the other two species.Another hypothesis is that the demographic histories of dung beetles in Mongolia were influenced more by late Quaternary climate change than by domestication events.Several previous studies have reported the influence of the late Quaternary climate oscillation, usually the last glacial maximum, on the population demography of species in Mongolia 35,36 .
The SDMs used in this study suggest that the distributions of dung beetles inhabiting a single biome type will contract and shift substantially in response to climate change.In contrast, the distributions of species inhabiting multiple biomes will expand, even if affected by range shifts.These trends are consistent with the idea that generalists and specialists respond differently to rapid environmental change.Species with broad tolerance exhibit range expansion, whereas species with restricted tolerance undergo range contraction 37,38 .
We revealed that the number of haplotypes would decrease following the extinction of local populations of all three species (Table 3).Under the most severe climate scenario (SSP585), C. erraticus and G. mopsus were predicted to experience haplotype losses of 54.5 and 68.9%, respectively, by 2070.Meanwhile, C. eumenes was predicted to have lost all extant haplotypes from the 10 sites investigated in this study.A reduction in the absolute number of haplotypes within the overall gene pool of a species is likely to reduce its evolutionary potential and resilience.However, despite reduced haplotype richness, the dung beetles modeled in this study were predicted to maintain genetic diversity (Table 3).These unexpected results may be attributed to the low mtDNA divergence of the species in both geographic configurations and lineages.In the phylogenetic trees of these species, unique haplotypes were scattered without geographic configuration or deep phylogenetic divergence.Therefore, haplotype loss due to local extinction would be significant under climate change scenarios, while low genetic divergence and a high singleton frequency could buffer the effects of climate change on genetic diversity and phylogeographic structure.Similar predictions have been reported by previous studies 19,39 .
Our predictions of future genetic diversity and distribution should be interpreted cautiously.First, the study could not address the extent to which dispersal might mitigate genetic diversity loss, and it was conservatively assumed that dispersal was restricted.In this study, local populations were assumed to disappear as habitats became climatically unsuitable, even when they were close to areas projected to persist or become climatically suitable.This assumption can enhance the prediction of potential genetic diversity loss, thereby contributing to effective conservation planning.However, it can reduce the prediction accuracy.For instance, the range of each species could shift due to enhanced fitness in locations at the leading edge and reduced fitness in locations at the trailing edge.This process tends to lead to a gradient in genetic diversity that increases toward the leading edge 32,40 .Second, the binary decision of presence or absence in the model neglected phenotypic plasticity and local adaptation.Several studies have investigated the local adaptation and phenotypic plasticity of dung beetles in terms of climatic differences.For example, Lim et al. (2020) reported variations in the morphology (body shape and size) of G. mopsus from steppe and desert in Mongolia 41 .Furthermore, Cuesta et al. 42 reported the phenological plasticity of Iberian dung beetles, and Macagno et al. (2018) reported the behavioral plasticity of the model species Onthophagus taurus, which was observed to mitigate the negative effects of temperature increases 43 .This study indicates that C. erraticus and G. mopsus inhabit climatically unsuitable localities (C.erraticus: CER1822 and CER1005; G. mopsus: GM1403, GM1510, and GM1630) with relatively high genetic diversity (except for CER1005, which included only one individual), which might suggest that, in Mongolia, the response of dung beetles to climatic stress involves either local adaptation or phenotypic plasticity.
Researchers have recently attempted to incorporate evolutionary factors into SDMs 44,45 .For example, Benito-Garzόn et al. 46 described a newly developed SDM that can incorporate fitness-related traits to estimate the sensitivity of species to environmental change.Abreu-Jardim et al. 39 forecasted the effects of climate change on the phylogeographic diversity of species.Nevertheless, such efforts remain in the early stages of development and have been applied to a few well-documented species.Determining the actual levels of genetic parameters required for responding to environmental change is challenging.Nonetheless, it's essential to conserve a wide range of genotypes across both phylogeny and geography to retain a meaningful subset of this diversity 21 .Our results highlight that the spatial arrangement of haplotypes and lineages, reflecting the distinct phylogeographic histories of each species, offers valuable insights into their resistance against potential genetic diversity loss owing to climate change.
As a result of climate change, the distribution of steppe-inhabiting dung beetles was expected to shift toward the mountainous steppe.In contrast, the distributions of desert-inhabiting dung beetles were expected to shift toward the central region and to expand northward toward the Gobi Desert (Fig. 5).However, contrary to our expectations, the desert-inhabiting species C. eumenes experienced the most severe range contraction.This suggests that the modeled climate change may have exceeded the thermal and drought tolerance of dung beetles in the Gobi Desert.Widespread dung beetles were predicted to lose their current habitat in the southern desert region and to expand their distribution into the northern-central region, which was predicted to develop into a more suitable habitat (i.e., semi-desert habitat).Therefore, it is likely that the central region of Mongolia, which is an intermediate region between desert and steppe, will become a key region for dung beetle diversity conservation over the next 50 years (Fig. 5).National parks in central Mongolia, such as Hustai, Gorkhi-Terelj, and Khugu Tarna, serve as climate refuges and corridors for Mongolian dung beetles.These national parks will likely become the leading edge of the range-shifting populations of C. eumenes and the trailing edge (or refuge) and center (or corridor) populations of C. erraticus and G. mopsus, respectively.Therefore, conservation management in the central region of Mongolia is crucial for conserving biodiversity of Mongolian dung beetles.
In conclusion, we provide novel insights into the potential range shifts and genetic diversity losses of dung beetles from the Gobi Desert and Mongolian Steppe in response to climate change.Contrary to general expectations, desert species appear to be the most vulnerable to climate change due to the extensive degradation in the Gobi Desert.Both steppe and widespread species are likely to experience high rates of occupancy turnover and genetic loss, which could reshuffle the dung beetles community.Climate change-induced desertification should be considered an emerging threat to dung beetles in drylands worldwide.

Genetic diversity and haplotype network analyses
The number of haplotypes (Nh), haplotype diversity (Hd), and nucleotide diversity (π) of each population were estimated for each species using DnaSP v.6.12.03 51 .Allelic richness (Ar) was calculated using the rarefaction method in CONTRIB v.1.025 52, that corrects for unequal sample sizes among populations.Haplotype networks for each species were constructed using minimum-spanning tree algorithms in HapStar v.0.7 53 and optimized using PopART v.1.7 54 .Only populations from which at least five individuals (N ≥ 5) were sampled were included in subsequent genetic analyses, except for demographic history and intraspecific phylogenetic analyses.

Population structure and demographic history analyses
To infer the genetic differentiation and spatial genetic structure of the three species, ARLEQUIN v.3.5 55 was used to calculate fixation index (F ST ) values for pairwise population comparisons and to conduct hierarchical analysis of molecular variance (AMOVA, 1000 permutations).Each population of the three species was assigned to a geographic or climatic region, depending on the species' distribution pattern.).The total molecular variance was partitioned among the groups (Φ ct = ′′inter-group′′ genetic variation), within the groups (Φ sc = ′′intra-group′′ genetic variation), and within populations regardless of the grouping (Φ st = ′′inter-population′′ genetic variation).
To infer demographic history, the past demographic expansion of the population of the three species was tested Tajima's D 56 , and Fu's Fs 57 neutrality tests with 1000 random samples in ARLEQUIN v3.5.Mismatch distribution (distribution of pairwise sequence difference) analysis was performed in DnaSP v.6.12.03 and sum of squared deviations (SSD) and raggedness index values were calculated to assess the goodness of fit between the observed and expected distributions of demographic expansion model.

Intraspecific phylogeny
Phylogenetic trees of the three species were reconstructed using Bayesian inference (BI) and Maximum likelihood (ML).BI trees were reconstructed using Markov Chain Monte Carlo (MCMC) methods, assuming a constant coalescent size for the tree prior, in BEAST v.10.4 58 .The best-fit evolutionary substitution models were selected using the Akaike information criterion (AICc) in jModelTest v.2.1.7 59 as follows: TrN + I for C. erraticus and C. eumenes; GTR + G + I for G. mopsus.Strict clock model was selected for three species based on calculation of the Bayes factor 58,60 .A COI substitution rate of 0.0075 (substitutions/site/million years) was based on previous phylogenetic study on dung beetles 61 .The MCMC run was performed using 10-100 million generations for each species, with sampling every 1000 generations, and the first 10 or 20% of generations were discarded as burn-in (Table S5).The run was examined using Tracer v.1.6 62to confirm that stationarity and convergence had been reached (effective sample size [ESS] > 200).The Maximum Clade Credibility (MCC) tree was evaluated using TreeAnnotator v.10.4 58 .ML analyses were performed using IQ-TREE 63 with 1000 replicates in ultrafast bootstrap 64 .Intraspecific coalescent trees were visualized using FigTree v.1.4.4 65 .

Species distribution modeling
From 2009 to 2022, field investigations were conducted at 164 sites distributed across Mongolia, ranging from 50° N to 43° N and 91° E to 112° E. Many SDM studies rely on distribution records from open-access databases, such as the Global Biodiversity Information Facility (GBIF).However, such databases often lack data for neglected regions and species, such as dung beetles in the Gobi Desert and Mongolian Steppe.To address this challenge, this study used data collected during our extensive field investigations conducted over the past 14 years (Table S6).Duplicate occurrence records were eliminated and rarefied to 5 km, resulting in 56, 28, and 45 occurrence records for C. erraticus, C. eumenes, and G. mopsus, respectively.
Bioclimatic variables for modeling (Table S7) were obtained from WorldClim 66 at a spatial resolution of 30 arcsec (ca. 1 km 2 ) using ArcGIS Desktop 10.8 67 .For future climate projections to 2050 and 2070, three different Shared Socioeconomic Pathways (SSPs), including two moderate stabilization scenarios (SSP245 and SSP370) and one high baseline emission scenario (SSP585), were derived from the CMIP6 global climate model (BCC-CSM) 68 .The limitation of food resources is a crucial factor in determining species' distribution 69 .However, in Mongolia, the tremendous number of livestock widely distributed through nomadic herding provides sufficient food resources across the landscape (National Statistics Office of Mongolia; https:// www.en.nso.mn/, accessed on April 08, 2022) 28,41 .
The current and future distributions of the three dung beetles were inferred using maximum entropy modeling 71 .The maximum entropy model (MaxEnt) is one of the most popular SDMs and uses presence-only data to estimate species distribution and habitat suitability 72 .The ′′ENMeval′′ R package 73 was used to select the best models for each species based on the optimization of two parameters: regularization multiplier (RM) and feature combination (FC).SDMs were generated using tenfold cross-validation with 10,000 background points.The models were optimized using RM values ranging from 0.5 to 4.0 (in increments of 0.5) and six feature combinations (L, LQ, H, LQH, LQHP, and LQHPT), which resulted in a total of 48 parameter combinations of SDMs for each species.The k-fold cross-validation approach was used instead of the delete-one-jackknife approach, which performs better for small sample sizes 74,75 .The optimized SDMs were evaluated using the difference between the training and testing AUCs (AUC diff ) and 10% training omission rate (OR 10 ), which indicate the degree of overfitting, and the Akaike information criterion corrected for small sample sizes (AICc), which reflects model complexity 14,74,76 .The optimal parameter combination for each species was selected based on the lowest delta.AICc value.To maintain a balance between model complexity and predictive accuracy, the values of OR 10 , AUC diff , and AUC test were considered collectively 14,77 .

Predicting distribution and genetic diversity under climate change scenarios
To simulate the future distributions and genetic diversities of the three dung beetle species under three climatic scenarios (SSP245, SSP370, and SSP585) at two time points (2050 and 2070), the predicted suitability across Mongolia was transformed into binary maps (presence/absence) based on the 10% training presence logistic threshold on the established final SDMs for each species.Changes in distribution were calculated as the difference in range size between current distributions and those predicted at each time point under each climatic scenario.We assumed that only populations occurring in areas with suitability greater than the present threshold would persist and, thus, contribute to the gene pool of subsequent generations 39 .The number of haplotypes, haplotype diversity, and nucleotide diversity were only calculated for the populations remaining at each time point under each climatic scenario.The loss of haplotypes under different scenarios was visualized in the phylogenetic tree for each species.

Figure 2 .
Figure 2. Distribution and environmental suitability predictions of Colobopterus erraticus.(A) Current species distribution model with the occurrence records (white dots), and sites for population genetics (black dots); (B-D) Predicted future distributions in 2070 under the SSP245, SSP370, and SSP585 scenarios, respectively.The map of distribution model was modified using ArcGIS 10.8 (www.esri.com) 67 .

Figure 3 .
Figure 3. Distribution and environmental suitability predictions of Cheironitis eumenes.(A) Current species distribution model with the occurrence records (white dots), and sites for population genetics (black dots); (B-D) Predicted future distributions in 2070 under the SSP245, SSP370, and SSP585 scenarios, respectively.The map of distribution model was modified using ArcGIS 10.8 (www.esri.com) 67 .

Figure 5 .
Figure 5. Future distribution of three model dung beetles in 2070 under the SSP585 scenario with schematic illustration of climate change and species response in the Gobi Desert and Mongolian Steppe.The map of distribution model was modified using ArcGIS 10.8 (www.esri.com) 67 .

Table 2 .
Analysis of molecular variance (AMOVA) for sampled populations of three dung beetle species in Mongolia.Significant values (P < 0.05) are indicated in bold.

Table 3 .
Predicted future genetic diversity and distribution for sampled populations of three dung beetle species in Mongolia under current conditions and future climate change scenarios.SpeciesScenarios/